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How the PHMC algorithm samples configuration space * 
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We show that in practical simulations of lattice QCD with two dynamical light fermion species the PHMC 
algorithm samples configuration space differently from the commonly used HMC algorithm. 



1. Introduction 

Simulations of lattice QCD with light dynam- 
ical fermions still pose severe problems, concern- 
ing both the computational cost and ergodicity 
properties. The PHMC algorithm Q represents 
an attempt of improving on both points over the 
commonly used HMC algorithm . 

The main idea [§ of the PHMC algorithm is 
to sample gauge configuration space -through a 
standard HMC method- by using an approxi- 
mate action where, in the fermion sector, the 
inverse of the squared Hermitian Dirac operator 
{Q^ oc M'^M) is replaced by a suitable polynomial 
approximant, Pn,e{Q^)- The use of a well con- 
trolled polynomial approximation makes also pos- 
sible to correct for it through an efficient reweigh- 
ing technique ||l[], leading to exact (reweighted) 
sample averages for lattice QCD with n/ = 2 de- 
generate fermion species: 



(O) - {W)p],,,c{OW)pHMC 



(1) 



where O stands for any observable and is a 
noisy estimate of det[Q^ Pn,e{Q^)]- 

According to the Chebyshev approximation 
method, the polynomial in s having degree n, 
Pn,e{s), approximates the function 1/s, with s > 
0, with a relative error that is bounded by (5 ~ 
2 exp(— 2-y/en) in the range|^ e < s < 1 and that 
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'^The operator is assumed to be normalized so that its 
highest eigenvalues is always smaller than 1. 



quickly increases as s gets smaller than e. As 
a consequence, in the molecular dynamics (MD) 
update, the role of the lowest eigenvalues of Q^, 
denoted here by Amin, is taken by e, which can 
be chosen -in practice- about 2(Aniin). This leads 
to a computational cost per MD trajectory cor- 
respondingly smaller -by about a factor 2- than 
in the case of the HMC algorithm. Moreover, 
the approximate action used in the MD evolution 
makes gauge configurations with very low (com- 
pared to e) eigenvalues of to be generated with 
much higher probability than in the HMC algo- 
rithm. These gauge configurations are expected 
to be important for the sample average of many 
fermion observables and for associated changes of 
topological sectors. We then expect the PHMC 
algorithm to show -in critical situations- ergodic- 
ity properties that are different from the ones of 
the HMC algorithm. 

2. Main results of PHMC tests 

The PHMC algorithm has been implemented 
on APE computers and carefully studied 
particular care has been devoted to the role of 
reweighing and to the tuning of n and e, for 
which a practical procedure was sugg ested 1). 
All test studies were performed on lattices with 
Schrodinger functional (SF) boundary conditions 
1^, which enabled us to work at vanishing quark 
mass. We have used Wilson fermions -in both 
the standard and the 0(a) improved versions- 
with even-odd preconditioning and a Sexton- 
Weingarten integration scheme for the MD evo- 
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Table 1 

The average condition number, (fc), of the pre- 
conditioned fermion matrix and the costs, Cq0, of 
a MD trajectory in performance tests at k k^. 
More details in refs. [1] (test a) and [5] (tests b,c). 



Lattice 












a: 8* 


5.6 





720 


7398 


3974 


b: 8^ • 16 


6.8 


1.4251 


760 


7750 


5956 


c: 8^ • 16 


5.4 


1.7275 


1500 


19734 


11450 



lution. Within the statistical uncertainties, we 
have always found consistent results for the mean 
values of several pure gauge and fermion observ- 
ables obtained from the HMC and the PHMC al- 
gorithm. 

Performance test^ against the HMC algorithm 
have been performed at vanishing quark mass on 
lattices of small and intermediate physical size, 
with spatial length never larger than 1 fm. This 
physical situation -although very different from 
the ones where most unquenched simulations are 
performed- is of practical interest for non pertur- 
bative renormalization studies. We summarize 
in table ^ the main results of our performance 
tests: the computational costs Cq^'^ and Cq™'^ 
are given in units of fermion matrix [Q) times 
pseudo-fermion vector ((^) multiplications and re- 
fer to a full trajectory, accounting in the PHMC 
case also for the reweighing procedure. However 
only for tests a and b comparing these costs corre- 
sponds to a comparison of the actual costs to gen- 
erate an independent gauge configuration, since 
for almost all the considered observables compat- 
ible errors (within 0(15%) relative uncertainties) 
are obtained from the two algorithms when the 
same statistics is employed. In the case of test c 
the accumulated statistics is not enough to make 
any definite statement on statistical errors and 
corresponding uncertainties. Moreover a direct 
performance comparison is made problematic by 
the significant differences in sampling the config- 
uration space that we are going to discuss. 

As shown in fig. 3 of ref. Q and fig. 4 of ref. ||], 
the PHMC sample in the case of test c includes 
a significant fraction of configurations carrying 
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We omit to discuss here many tests on 4 lattices 



Figure 1. An example of Monte Carlo time evo- 
lution for Amin and the naive topological charge 
after cooling, Qtopo- 

one (or few) isolated modes of the fermion ma- 
trix lying some orders of magnitude below the 
average value of Amin- Moreover, different values 
of the naive topological charge are measured after 
500 cooHng iterations, suggesting that changes of 
topological sectors do occur - even at zero quark 
mass and in a space volume less than 1 fm^. A 
typical example is shown in fig. |l|. We recall that 
no index theorem has to hold in this case and 
refer to |^ for further considerations. No excep- 
tionally small eigenvalues were observed in the 
corresponding HMC simulation . 

We present in fig. |^ the Monte Carlo time evolu- 
tion of the four-fermions Green function fA{T/2) 
(see e.g. eq.(14) of ref. for its definition), 
taking two typical Monte Carlo history segments 
from our PHMC and HMC data. In the PHMC 
case we see that peaks for /^(T/2) occur in coin- 
cidence with very small values of Amin- The con- 
tributions to QCD sample averages, eq.(^, com- 
ing from these "exceptional" configurations, are 
made of "normal" size by the corresponding small 
values of W . Many examples of this behaviour, 
with even larger spikes in the values of /^(T/2) 
and Amin, were seen in the PHMC data. A sim- 
ilar behaviour has been also observed in numer- 
ical studies of Supersymmetry ||]. On the other 
hand, the HMC algorithm seems to generate with 
very low probability these exceptional configura- 
tions, from which a relevant finite contribution to 
QCD sample averages for many fermion observ- 
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Figure 2. A segment of the Monte Carlo history of /yi(T/2) as obtained from the PHMC and the HMC 
algorithms. In the PHMC case we also show the corresponding history of Amin(Q^) and W. 



ables may in principle come. 

3. Conclusions 

We have shown that the PHMC algorithm, as 
expected, samples the gauge configuration space 
differently from the HMC algorithm, while be- 
ing at least competitive with it from the per- 
formance point of view. If gauge configurations 
carrying exceptionally small eigenvalues of the 
fermion matrix are important for some observ- 
ables, the PHMC algorithm should be largely su- 
perior, owing to its ability^ of sampling and prop- 
erly treating these configurations. We think that 
consequences of the different ergodicity and per- 
formance properties of the HMC and the PHMC 
algorithm deserve to be further studied with more 
statistics and on larger volumes. 

This work is part of the ALPHA collaboration 
research programme. We thank DESY for allo- 



^ This remains true even in presence of exact fermion zero 
modes, as discussed in M. 



eating computer time to this project. 
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